home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / RNDGEN10 / RANGEN.PAS < prev    next >
Pascal/Delphi Source File  |  1993-12-06  |  12KB  |  380 lines
  1. {$A-,B-,D+,E+,F-,I+,L+,N+,O-,R+,S+,V+}
  2. UNIT Rangen;
  3.  
  4. {==============================================================================
  5.  
  6. MODULE:         RANDOMGENERATOR (RANGEN.PAS)
  7.  
  8. ===============================================================================
  9.  
  10. VERSION:        V1.5            Turbo-Pascal 6.0, 7.0
  11.  
  12. COPYRIGHT:      (c) 1987-1993
  13.                 Dr. W. Gross, Keplerstrasse 51
  14.                 D-69120 Heidelberg, FRG
  15.                 gross@aecds.exchi.uni-heidelberg.de
  16.  
  17. CREATED:        10-NOV-87       Wolfgang Gross          on VAX/VMS, PASCAL-2
  18.  
  19. LAST UPDATE:    02-MAY-88       Wolfgang Gross          round off error in
  20.                                                         NORRANDOM
  21.                 21-JUL-88       Wolfgang Gross          round off error in
  22.                                                         UNIRANDOM
  23.                 15-MAR-89       Wolfgang Gross          NorRan, RanSphere
  24.  
  25. TESTED:         10-FEB-90       Wolfgang Gross
  26.  
  27. -------------------------------------------------------------------------------
  28. TITLE:          RANDOM NUMBER GENERATES
  29. -------------------------------------------------------------------------------
  30.  
  31. DESCRIPTION:    Random number generators for uniform, normal, exponential,
  32.                 lognormal, poisson and uniform sphere distribution.
  33.  
  34.                 If the TURBO-PASCAL random generator RANDOM is used
  35.                 one must set the variable SYSTEM.RandSeed before calling
  36.                 RANDOM and save it after the call using the given
  37.                 seed variable. This allows independent streams of
  38.                 random numbers. See comment in UNIRANDOM.
  39.  
  40.                 If TEST8086>2 we can use 386 instructions. UNIRANDOM
  41.                 will then use its own method to generate the next
  42.                 univariate pseudo random number as described below.
  43.                 The method is approved by many independent tests
  44.                 that have been published in the literature. It has
  45.                 good Bayer coefficents which means it can be used for
  46.                 generating multi-dimensional random vectors.
  47.                 One might also utilize a floating point processor
  48.                 or switch to IEEE reals (type single, double) etc.
  49.                 Everybody may elaborate on this him/herself.
  50.  
  51.                 If someone needs to modify and recompile this unit
  52.                 and doesn't have TASM, the complete ASM code is
  53.                 provided for compilation by the internal assembler.
  54.                 Some special arrangements are necessary to cope
  55.                 with 386 instructions. Define HAVETASM if you want
  56.                 to link MD2P31.OBJ. If this symbol is not defined
  57.                 the inline assembler version is used.
  58.  
  59.  
  60.  
  61.  
  62.                 This unit has evolved from a similar module worked
  63.                 out under PASCAL-2 in a VAX/VMS environment (including
  64.                 VAX assembler for MD2P31). The module can be found
  65.                 with anon ftp on aecds.exchi.uni-heidelberg.de.
  66.  
  67.  
  68. CALLED BY:      <module name>   <title>
  69.  
  70.  
  71. EXPORTS:        UNIRANDOM       uniformly distributed random numbers
  72.  
  73.                 NORRANDOM       normally distributed random numbers
  74.  
  75.                 NORRANDOM2      normally distributed random numbers,
  76.                                 polar method
  77.  
  78.                 NORRANDOM3      normally distributed random numbers,
  79.                                 other mean/std-dev values than (0,1)
  80.  
  81.                 EXPRANDOM       exponentially distributed random numbers
  82.  
  83.                 LOGNORRANDOM    logarithmic normal distribution
  84.  
  85.                 POIRANDOM       Poisson distribution
  86.  
  87.                 RanSphere       uniformly distri. random numbers on sphere
  88.  
  89. INPUT:          seed            random number seed
  90.  
  91.                 mean            mean value for distribution
  92.  
  93.                 std             standard deviation
  94.  
  95. OUTPUT:         function value = random number
  96.  
  97. }
  98.  
  99.  
  100. {=============================================================================}
  101.  
  102. { Exports: }
  103. Interface
  104.  
  105.  
  106. FUNCTION  UNIRANDOM  ( VAR seed : longint ) : real;
  107.  
  108. FUNCTION  NORRANDOM  ( VAR seed : longint ) : real;
  109.  
  110. FUNCTION  NORRANDOM2 ( VAR seed : longint;
  111.                        VAR v2   : real )    : real;
  112.  
  113. FUNCTION  NORRANDOM3 ( mean, std: real;
  114.                        VAR seed : longint ) : real;
  115.  
  116. FUNCTION  EXPRANDOM  ( mean     : real;
  117.                        VAR seed : longint ) : real;
  118.  
  119. FUNCTION  LOGNORRANDOM ( mean, std : real;
  120.                        VAR seed : longint ) : real;
  121.  
  122. FUNCTION  POIRANDOM  ( lambda   : real;
  123.                        VAR seed : longint ) : longint;
  124.  
  125. PROCEDURE RanSphere  ( VAR seed : longint;
  126.                        VAR x,y,z: real );
  127.  
  128.  
  129. {===========================================================================}
  130.  
  131. Implementation
  132.  
  133.  
  134. {DEFINE HAVETASM} {put a $ in front of DEFINE to link MD2P31.OBJ}
  135.  
  136. {$IFDEF HAVETASM}
  137.  
  138.    PROCEDURE MD2P31 (A,B : longint;
  139.                      VAR Q,R : longint); far; external;
  140.  
  141.  
  142.    {$L MD2P31.OBJ}
  143.  
  144. {$ELSE}
  145.  
  146.    PROCEDURE MD2P31 (A,B : longint;
  147.                      VAR Q,R : longint); far; assembler;
  148.      {rather tricky: we fake TP into generating 386 intructions
  149.       Calculate product a*b, represent product as q*2^31+r, 0<=r<2^31}
  150.      asm
  151.        push  es
  152.        push  di
  153.  
  154.        db    $66       { 386 prefix for dword operation}
  155.        push  ax        { = push eax }
  156.  
  157.        db    $66
  158.        push  bx        { = push ebx }
  159.  
  160.        db    $66
  161.        push  dx        { = push edx }
  162.  
  163.        db    $66
  164.        mov   ax,word ptr ss:[bp+18]  { mov   eax,a }
  165.  
  166.        db    $66
  167.        mov   bx,word ptr ss:[bp+14]  { mov   ebx,b }
  168.  
  169.        db    $66
  170.        imul  bx        { imul ebx   ; a*b in edx:eax  (= 8 bytes !) }
  171.  
  172.        db    $66
  173.        mov   bx,ax     { mov ebx,eax }
  174.  
  175.        db    $66,$0f,$a4,$da,$01     { SHLD edx,ebx,1      ; edx contains q }
  176.  
  177.        db    $66,$25,$ff,$ff,$ff,$7f { and  eax,07fffffffh ; eax contains r }
  178.  
  179.        les   di,r
  180.        db    $66
  181.        stosw           { stosd   ; mov eax to r }
  182.  
  183.        db    $66
  184.        mov   ax,dx     { mov eax,edx }
  185.  
  186.        les   di,q
  187.        db    $66
  188.        stosw           { stosd   ; mov eax to q }
  189.  
  190.  
  191.        db   $66
  192.        pop  dx         { = pop edx }
  193.  
  194.        db   $66
  195.        pop  bx         { = pop ebx }
  196.  
  197.        db   $66
  198.        pop  ax         { = pop eax }
  199.  
  200.        pop  es
  201.        pop  di
  202.      end;{ PROC MD2P31 }
  203.  
  204.  
  205. {$ENDIF}
  206.  
  207.  
  208. FUNCTION UNIRANDOM;
  209.   { Cf. Payne,W.H., Rabung, J.R., Bogyo, T.P. :
  210.         "Coding the Lehmer pseudo-random number generator.
  211.     Comm. ACM 12, 85-86 (1969)                         }
  212.  
  213.   CONST MODULUS : longint = 2147483647;  { = 2^31-1 }
  214.         FACTOR  : longint = 397204094;   { primitive root, e.g. used by SAS }
  215.  
  216.   VAR   Q,R,S : longint;
  217.  
  218.   BEGIN { UNIRANDOM }
  219.  
  220.     IF TEST8086>1 THEN
  221.       BEGIN
  222.         { Priniple
  223.          seed := ( seed*FACTOR )  MOD  MODULUS;
  224.           but can not use it due to overflow }
  225.  
  226.         MD2P31 ( seed, factor, Q, R );
  227.         S := modulus - R;
  228.         IF S > Q
  229.           THEN seed := Q + R
  230.           ELSE seed := Q - S;
  231.  
  232.         { Single precision version. For more details on the division cf. to
  233.           Fishman, G.S., Moore, L.R. :
  234.           A statistical evaluation of multiplicative congruential random
  235.           number generators with modulus 2^31-1 .
  236.           Jour. Amer. Stat. Ass. 77, 129-136 *1982)                         }
  237.  
  238.         UNIRANDOM := seed/MODULUS;
  239.       END
  240.      ELSE
  241.       {see comment at the beginning of unit, RandSeed from unit SYSTEM}
  242.        BEGIN
  243.          RandSeed := seed; Unirandom := random; seed := RandSeed;
  244.        END;
  245.  
  246.  
  247.   END; { FUNC UNIRANDOM }
  248.  
  249.  
  250.             {--------------------------------------------}
  251.  
  252.  
  253. FUNCTION  NORRANDOM  ( VAR seed : longint  )  : real;
  254.   CONST TwoPi = 2*Pi;
  255.   VAR     s, r1, r2, z1, z2 : real;
  256.   BEGIN { NORRANDOM }
  257.     REPEAT r1 := UniRandom ( seed ) UNTIL r1>0;
  258.     r2 := UniRandom ( seed );
  259.     s := -2*ln(r1);
  260.     IF s<0 THEN s := 0;         { round off error could give negative value s }
  261.     NORRANDOM := SQRT ( s ) * cos ( TwoPi*r2 );
  262.   END; { FUNC NORRANDOM }
  263.  
  264.             {--------------------------------------------}
  265.  
  266. FUNCTION  NorRandom2 ( VAR seed : longint;
  267.                        VAR v2   : real ) : real;
  268.  { Polar method due to Box, Mueller and Marsaglia,
  269.    cf. D.E. Knuth, The art of computer porgramming, vol 2, 117
  270.    first call: v2 must be > 1000 }
  271.  
  272.   VAR     s, t, v1 : real;
  273.  
  274.   BEGIN { NORRANDOM2 }
  275.     IF v2<1000
  276.       THEN
  277.         BEGIN
  278.           NORRANDOM2 := v2; v2 := 2000
  279.         END
  280.       ELSE
  281.         BEGIN
  282.           REPEAT
  283.             v1 := UniRandom ( seed ); v1 := v1+v1-1;
  284.             v2 := UniRandom ( seed ); v2 := v2+v2-1;
  285.             s := Sqr (v1) + sqr(v2);
  286.           UNTIL (s<1) AND (s>0);
  287.           T := Sqrt ( -2*Ln(s)/s );
  288.           NORRANDOM2 := V1*T;
  289.           V2 := V2*T;
  290.         END;
  291.   END; { FUNC NORRANDOM2 }
  292.  
  293.  
  294.             {--------------------------------------------}
  295.  
  296. FUNCTION  NORRANDOM3 ( mean, std : real;
  297.                        VAR seed : longint ) : real;
  298.   BEGIN { NORRANDOM3 }
  299.     NORRANDOM3 := std*NORRANDOM (seed) + mean;
  300.   END; { FUNC NORRANDOM3 }
  301.  
  302.             {--------------------------------------------}
  303.  
  304. FUNCTION  EXPRANDOM  ( mean     : real;
  305.                        VAR seed : longint ) : real;
  306.   BEGIN { EXPRANDOM }
  307.      EXPRANDOM := ( -ln(UNIRANDOM(seed))*mean );
  308.   END; { FUNC EXPRANDOM }
  309.  
  310.             {--------------------------------------------}
  311.  
  312. FUNCTION  LOGNORRANDOM ( mean, std : real;
  313.                          VAR seed  : longint ) : real;
  314.   VAR     m2, s2, sum, mu, sigma : real;
  315.   BEGIN { LOGNORRANDOM }
  316.     m2 := SQR(mean); s2 := SQR(std); sum := m2+s2;
  317.     mu := 0.5*ln(SQR(m2)/sum);
  318.     sigma := SQRT (ln(sum/m2));
  319.     LOGNORRANDOM := exp (sigma * NORRANDOM (seed) + mu);
  320.   END; { FUNC LOGNORRANDOM }
  321.  
  322.             {--------------------------------------------}
  323.  
  324. FUNCTION  POIRANDOM  (  lambda   : real;
  325.                         VAR seed : longint )  : longint;
  326.  
  327.  
  328.   VAR     i : longint;
  329.           p,q,r,z : real;
  330.  
  331.   BEGIN { POIRANDOM }
  332.     { cf. H.E. Schaffer, Generator of random numbers satisfying the
  333.           Poisson distribution, Comm. ACM 13, 49 (1970)
  334.           G.S. Fishman, Sampling from the Poisson distribution on
  335.           a computer, Computing 17, 147-156 (1976) }
  336.     IF lambda<50
  337.       THEN
  338.         BEGIN
  339.           z := exp ( -lambda ); p := 1; i := -1;
  340.           REPEAT
  341.             i := i+1;
  342.             r := UNIRANDOM ( seed );
  343.             p := p*r;
  344.           UNTIL p<=z;
  345.           POIRANDOM := i;
  346.         END
  347.       ELSE
  348.         BEGIN
  349.           i := Round ( NORRANDOM3 ( lambda, SQRT(lambda), seed ) );
  350.           IF i<0 THEN i :=0;
  351.           POIRANDOM := i;
  352.         END;
  353.  
  354.   END; { PROC POIRANDOM }
  355.  
  356.             {--------------------------------------------}
  357.  
  358. PROCEDURE RanSphere ( VAR seed : longint;
  359.                       VAR x,y,z : real );
  360.  { R.E. Knop, Random Vectors Uniform in Solid Angle,
  361.    CACM 13 (1970), 326 }
  362.  
  363.   VAR     s, t, v1, v2 : real;
  364.  
  365.   BEGIN { RanSphere }
  366.  
  367.     REPEAT
  368.       v1 := UniRandom ( seed ); v1 := v1+v1-1;
  369.       v2 := UniRandom ( seed ); v2 := v2+v2-1;
  370.       s := Sqr (v1) + sqr(v2);
  371.     UNTIL (s<1) AND (s>0);
  372.  
  373.     T := 2*Sqrt ( 1-S );
  374.     x := T*V1; y := T*V2; z := S+S-1;
  375.  
  376.   END; { PROC RanSphere }
  377.  
  378.  
  379. end. {UNIT RANGEN}
  380.